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ABSTRACT 

This  report  describes  a  vibration  energy  harvesting  approach  that  uses  a  magnetoelectric  (ME) 
transducer  to  harvest  energy  from  bi-axial  vibrations.  The  approach  is  being  explored  as  a 
potential  means  of  powering  in  situ  structural  health  monitoring  systems  embedded  within 
aircraft  and  other  high  value  engineering  assets  that  experience  mechanical  vibration.  A  bi-axial 
oscillator  is  created  using  a  permanent-magnet/ ball-bearing  arrangement,  which  has  the  added 
benefit  of  permitting  a  relatively  compact  design.  The  magnet  produces  a  bi-axial  restoring  force 
on  the  bearing,  and  as  the  bearing  oscillates  it  steers  a  magnetic  field  through  a 
magnetostrictive/ piezoelectric  laminate  transducer  thereby  producing  an  oscillating  charge  that 
can  be  harvested.  A  simple  laboratory  demonstrator  of  a  bi-axial  ME  harvester  was  created  using 
a  Terfenol-D/lead  zirconate  titanate/Terfenol-D  transducer,  and  was  shown  to  produce  a  peak 
rms  power  of  121  pW  from  an  rms  acceleration  of  61  mG  at  9.8  Hz. 
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Bi-axial  Vibration  Energy  Harvesting 


Executive  Summary 

In-situ  Structural  Health  Monitoring  (SHM)  devices  allow  the  Australian  Defence  Force  to 
move  from  expensive  time-based  maintenance  approaches  for  ageing  platforms  to  cost- 
effective  condition-based  approaches.  For  air  platforms  the  installation  of  these  systems  is 
complicated  by  the  fact  that  the  majority  of  SHM  devices  need  to  be  fitted  on  internal 
aircraft  structure,  underneath  the  aircraft's  skin.  If  the  SHM  device  is  in  a  location  that  is 
difficult  to  access,  then  powering  the  device  may  be  problematic  because  traditional 
powering  methods  are  in  general  not  feasible.  For  example,  replacing  batteries  on  many 
SHM  devices  deployed  across  a  fleet  would  be  impractical,  and  accessing  an  on-board 
power  system  to  supply  SHM  devices  may  lead  to  flight  worthiness  and  certification 
issues.  To  address  this  powering  issue  the  Australian  Defence  Science  and  Technology 
Organisation  (DSTO)  is  investigating  the  possible  use  of  vibration  energy  harvesting 
(VEH).  Two  unresolved  scientific  issues  that  inhibit  the  use  of  VEH  on  aircraft  are:  (i)  the 
need  for  a  wide  operational  frequency  bandwidth  to  permit  harvesting  from  the 
frequency-rich  vibration  that  can  be  present  on  airframes,  and  (ii)  the  need  for  a  multi- 
axial  harvesting  approach,  since  aircraft  vibrations  are  typically  not  uni-axial.  Previous 
collaborative  work  between  the  DSTO  and  the  Active  Materials  Laboratory  at  UCLA 
addressed  the  first  issue  by  developing  the  vibro-impacting  energy  harvesting  approach 
which  produced  VEH  over  a  broader  operational  bandwidth  compared  with  many  other 
harvester  approaches,  including  harvesters  that  are  currently  commercially  available.  The 
second  fundamental  issue  with  most  VEH  approaches  (again  including  all  known 
commercial  vibration  energy  harvesters)  is  that  they  are  uni-directional,  and  hence  can 
only  harvest  vibrational  energy  from  host  accelerations  along  a  single  axis.  Therefore, 
while  a  considerable  amount  of  scientific  literature  exists  on  the  topic  of  VEH,  none  to  date 
reports  on  a  technique  to  effectively  harvest  from  bi-axial  host  accelerations.  This  report 
describes  a  bi-axial  approach  that  represents  a  significant  advancement  in  VEH, 
specifically  the  approach  increases  the  operational  directionality  from  single-axis  to 
360  degrees  in  a  plane.  Furthermore,  to  the  author's  knowledge  this  is  the  first  harvester 
design  that  uses  a  magnet/ bearing  cantilever  analogue  (replacing  the  cantilever  design 
used  by  many  harvesters  described  in  the  literature)  potentially  allowing  a  significant 
reduction  in  harvester  volume.  Finally,  to  the  authors'  knowledge  the  harvester  described 
in  this  report  is  the  first  that  uses  an  oscillating  ball-bearing  to  create  magnetic  flux 
steerage  through  a  magnetoelectric  laminate  transducer  to  generate  harvestable  electrical 
power.  This  report  will  describe  modelling  of  the  bi-axial  harvester,  and  will  also  report  on 
a  simple  laboratory  demonstrator  that  was  developed  as  a  proof  of  concept. 
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1.  Background 

The  Australian  Defence  Science  and  Technology  Organisation  (DSTO)  is  developing  a  variety 
of  in-situ  structural  health  monitoring  (SHM)  approaches  [1, 2, 3, 4, 5]  for  potential  use  in  high 
value  platforms  across  the  Australian  Defence  Force.  The  SLIM  systems  under  development 
could  be  employed  to:  (i)  continuously  monitor  airframe  loads  and  accelerations  during  flight, 

(ii)  detect  damage  and  damage  growth  and  other  structural  problems,  and  (iii)  provide  a  basis 
for  near-real-time  damage  assessment.  This  technology  could  also  potentially  permit  a  safe 
reduction  in  inspection  and  regular  maintenance  costs  and  therefore  reduce  aircraft  through- 
life  support  costs. 

DSTO  is  currently  investigating  the  various  components  of  the  generic  SLIM  concept  depicted 
in  Figure  1.  The  concept  involves  three  main  components,  being:  (i)  a  sensor  mounted  inside 
the  aircraft  at  a  difficult-to-access  location,  that  is  monitoring  in-flight  mechanical  loads  on  an 
airframe  [5],  (ii)  with  the  sensor  utilising  energy  that  is  parasitically  harvested  from  local 
airframe  vibrations  by  an  energy  harvester  [6,  7,  8,  9],  and  (iii)  when  the  aircraft  is  on  the 
ground  a  wireless  link,  the  acoustic  electric  feedthrough,  is  used  to  download  sensor  data  and 
to  simultaneously  provide  additional  energy  to  the  sensor  unit  [10, 11, 12, 13, 14]. 


•  Sensing  Points 


Harvester  Feedthrough 

Figure  1.  Schematic  of  a  wireless  Structural  Health  Monitoring  system  concept  with  sensing  unit, 
energy  harvester  and  wireless  data  and  power  transfer  capability 
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1.1  Magnetostriction,  piezoelectricity,  the  magnetoelectric  effect,  and 
magnetostrictive/piezoelectric  lamina  composites 

In  1847  Joule  observed  that  "the  elongation  is  proportional,  in  a  given  [ferromagnetic]  bar,  to 
the  square  of  the  magnetic  intensity"  [15],  the  behaviour  was  termed  the  'Joule  Effect'  and 
became  known  as  magnetostriction.  Somewhat  later  the  reciprocal  behaviour  (or  the  inverse 
magnetostrictive  effect)  was  discovered  by  Villari,  who  found  that  when  a  mechanical  stress  is 
applied  to  a  ferromagnetic  material  that  a  change  in  magnetisation  is  produced  [16].  In  the 
1970's  a  magnetostrictive  alloy  with  the  formula  TbxDyi_xFe2  (x  ~  0.3)  was  discovered  by 
researchers  at  the  U.S.  Naval  Ordinance  Laboratory.  The  alloy  was  named  Terfenol-D  after  the 
laboratory  (nol)  and  also  its  constitutive  elements,  terbium  (ter),  iron  (fe)  and  dysprosium  (D). 
The  alloy  exhibits  what  is  known  as  the  'giant  magnetostrictive  effect',  producing  strains 
hundreds  of  times  larger  that  those  of  classical  magnetostrictive  materials  such  as  iron  [17]. 

In  1880  the  Curie  brothers  demonstrated  the  direct  piezoelectric  effect,  finding  that  an  applied 
mechanical  stress  could  produce  electrical  charge  in  a  solid  [18].  The  reciprocal  behaviour  (or 
the  reverse  piezoelectric  effect)  was  mathematically  deduced  by  Lippman's  discussion  of 
fundamental  thermodynamic  conservation  laws  in  reference  [19].  Physicists  at  the  Tokyo 
Institute  of  Technology  investigated  various  piezoceramic  materials  [20],  developing  lead 
zirconate  titanate  (PZT)  in  around  1952.  Reference  [16]  describes  PZT  as  "the  most  common 
ferroelectric  compound  employed  for  smart  material  applications",  in  part  because  of  its  high 
electromechanical  coupling  and  extended  operating  temperature  range. 


(a) 


Polarization 
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— *  «— 
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Figure  2.  Magnetoelectric  lamina  composites  using  Terfenol-D  and  PZT  disks,  as  described  by  Ryu  et 
al  [25],  (Reproduced  with  permission,  Japanese  Society  of  Applied  Physics). 

The  linear  magnetoelectric  (ME)  effect  describes  the  induction  of  magnetisation  by  an  electric 
field,  or  the  reverse  effect  of  electric  polarisation  by  a  magnetic  field.  The  history  of  discovery 
of  the  ME  effect  is  well  described  in  reference  [21].  The  effect  was  termed  such  by  Debye  in 
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1926  [22]  after  earlier  work  by  Roentgen  [23]  and  Curie  [24],  In  2001  Ryu  et  al  [25]  developed  a 
magnetostrictive/ piezoelectric  lamina  composite  with  a  ME  voltage  coefficient  of  4.68  V  cm 1 
Oe1;  the  coefficient  was  36  times  greater  than  any  previously  reported  ME  particulate 
composite.  Figure  2  is  reproduced  from  Ryu's  article,  with  Figure  2a  depicting  the  behaviour 
of  a  Terfenol-D/PZT  laminate  under  an  applied  magnetic  field.  Reference  [25]  develops  an 
expression  for  the  (open  circuit)  output  voltage  generated  by  an  ME  transducer  subject  to  a 
stress  in  the  direction  orthogonal  (i.e.  the  1-direction,  or  y-axis)  to  the  direction  of  polarisation 
(i.e.  the  '3-direction',  or  z-axis), 

VoUT  =2x  g^l  x  tp  X  CT^ip  ,  (1) 

where  g3i  is  the  piezoelectric  voltage  constant  (Vm/N),  tp  is  the  thickness  of  the  PZT  layer  (m), 
and  <j^p  is  the  mechanical  stress  in  the  PZT  1-direction  (Pa). 

1.2  Vibration  energy  harvesting 

The  use  of  parasitically  harvested  energy  to  power  small  sensor  systems  has  become  viable 
over  the  past  decade  due  to  the  ever-decreasing  power  requirements  of  electronic  devices  [26]. 
This  has  resulted  in  considerable  scientific  interest  in  energy  harvesting  technologies  [27, 28]. 

In  particular,  there  has  been  significant  interest  in  the  area  of  vibration  energy  harvesting 
(VEFI)  [29],  particularly  for  single-degree-of-freedom  [30]  (SDOF)  piezoelectric  [31]  and 
magnetic  [32]  harvesters.  Although  there  are  a  number  of  commercially  available  SDOF 
harvesters  [27],  there  exist  a  number  of  continuing  challenges  that  prevent  widespread 
application  of  VEH  such  as  limited  operational  frequency  bandwidth  [9]  and  low  power 
density  [33] .  In  many  practical  applications  the  direction  of  the  incident  vibration  may  vary,  so 
another  challenge  is  the  uni-axial  nature  of  many  harvesters  [34]  particularly  cantilever-based 
piezoelectric  designs  [35]  and  SDOF  electromagnetic  designs  [36],  The  magnetoelectric  (ME) 
effect  produced  via  magnetostrictive/ piezoelectric  laminates  [37]  has  been  used  recently  in  a 
number  of  VEH  studies  [38].  In  particular,  a  recent  investigation  [39]  employed  a  cantilever 
arrangement  with  a  tip-mass  consisting  of  two  pairs  of  rare-earth  (NdFeB)  magnets,  that 
oscillates  around  a  fixed  Terfenol-D/PbTi03  laminate  transducer.  The  uni-axial  harvester 
examined  in  reference  [39],  with  a  single  laminate  transducer,  produced  -100  pW  from  a 
-33  Hz,  100  mG  excitation  (assumed  to  be  an  rms  acceleration). 

This  report  describes  an  inherently  compact  bi-axial  VEH  approach  that  has  been  developed 
by  DSTO  [40].  The  approach  uses  a  permanent-magnet/ ball-bearing  combination,  that 
operates  in  conjunction  with  a  Terfenol-D/PZT  laminate  transducer  (hereafter  called  the  'ME 
transducer').  On  the  4th  May  2011,  the  DSTO  secured  U.S.  provisional  patent  protection  on  the 
bi-axial  harvesting  approach  described  in  this  report  [41], 
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2.  Modelling 


This  section  examines  the  two  and  three-dimensional  models  developed  to  investigate  the  bi¬ 
axial  harvester  described  in  reference  [41],  which  describes  three  nominal  harvester  types. 
This  report  will  detail  one  of  the  harvester  types,  which  has  an  ME  transducer  located 
between  a  ball-bearing  (hereafter  called  the  'bearing')  and  a  magnet.  An  artist's  impression  of 
what  this  type  of  harvester  might  look  like  in  practice  is  shown  in  Figure  3. 
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Shown  in  (a)  is  an  exploded  view  of  the  hypothetical  harvester,  and  in  (b)  a  view  of  the 
assembled  harvester. 
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The  bearing  shown  in  Figure  3  is  subject  to  a  magnetic  restoring  force  in  the  x-y  plane, 
however  is  otherwise  free  to  move  across  the  surface  of  the  bearing  housing  in  response  to 
host  structure  oscillations.  Comsol  [42]  multi-physics  software  was  used  to  make  predictions 
about  the  restoring  force,  and  also  about  the  magnetic  field  distribution  as  the  bearing  location 
changes. 

2.1  Two  dimensional  modelling 

A  two-dimensional  model  was  created  based  on  the  arrangement  shown  in  Figure  4.  This 
arrangement  was  chosen  since  it  reflected  the  geometry  of  the  device  that  was  built  to 
demonstrate  the  principle  of  the  bi-axial  ME  harvester  (the  device  will  be  described  in  detail  in 
Section  3  of  this  report).  Note  that  the  central-line  of  the  harvester  is  defined  as  the  vertical 
line  (z-direction)  through  the  centre  of  the  magnet  (as  shown  in  Figure  4).  Rotational  effects 
were  ignored. 


Z 


y 


A 


ME  Transducer 


B 


Ball  bearing  (steel  CT-plated) 

Magnetostrictive  (Terfenol-D) 

Piezoceramic  (PZT-5A) 

Magnets  (NdFeB) 

Host  structure  (Fe) 

z-axis  central-line 


Figure  4.  Two-dimensional  schematic  diagram  of  the  bi-axial  energy  harvester,  showing  the 
magnetoelectric  (ME)  transducer  located  between  the  bearing  and  the  magnet.  Also  shown 
is  the  bearing  in  its  central  rest  position  (A)  and  the  magnetic  restoring  force  Fy,  acting  on 
the  displaced  bearing  (B).  The  bearing  has  a  diameter  of  25. 4  mm.  (Originally  published  in 
reference  [40]). 


The  model  was  analysed  using  Comsol,  allowing  predictions  of  the  static  magnetic  field 
distribution  to  be  made  (for  various  bearing  positions).  Of  particular  interest  is  the 
distribution  of  magnetic  flux  B  in  the  ME  transducer  since  this  determines  the  stress 
distribution  in  the  transducer,  and  hence  the  voltage  generated.  The  Comsol  application  mode 
‘perpendicular  induction  currents,  vector  potential'  was  used  to  calculate  the  static  magnetic  field 
distribution  for  various  bearing  locations.  In  particular,  the  calculated  magnetic  fields  were 
obtained  from  the  static  partial  differential  formulation  of  Ampere's  law  [43], 
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V  x(//01  V  x  A-M)-<jVV =Je , 


(2) 


where  M  is  the  in-plane  (y-z)  magnetisation  vector,  A  is  the  magnetic  vector  potential,  V  is  the 
electric  potential,  o  is  the  electric  conductivity,  go  is  the  permeability  of  a  vacuum,  and  Je  is 
the  vector  current  density. 


Two-dimensional  modelling  of  the  Terfenol-D  nonlinear  magnetostrictive  response  was 
carried  out  using  an  approach  adapted  from  reference  [44],  which  uses  the  phenomenological 
magnetostrictive  equations  [45, 46] .  The  deformation  due  to  magnetostriction  can  be  obtained 
using  the  following  equation. 


SI 

/ 


t.  f  1 A 

—X.  cos2  9  — 

2  \  3) 


(3) 


where  l  is  the  length  of  magnetostrictive  material,  Xs  is  a  constant  representing  the 

magnetostrictive  deformation  at  saturation,  and  0  is  the  angle  between  the  directions  of 
magnetisation  and  deformation.  Lee,  in  reference  [45],  states  that  equation  (3)  "seems  to  hold 
fairly  well  for  polycrystalline  samples  containing  no  preferred  orientation  of  crystal  axes". 
Polycrystalline  Terfenol-D  satisfies  this  criterion  since  it  is  nominally  isotropic.  The  directional 
term  can  be  rewritten  as  the  fraction  of  magnetisation  with  respect  to  saturation  magnetisation 
thus. 


-  =  -X 

l  2  5 


(4) 


where  I  is  the  magnetisation  and  I s  is  the  magnetisation  at  saturation.  The  1/3  term  takes  into 
account  the  magnetisation  required  to  align  the  individual  domains  within  the 
magnetostrictive  material,  and  for  this  model  it  is  ignored.  This  assumes  that  all  domains  are 
aligned,  and  furthermore  are  aligned  perpendicular  with  the  direction  of  magnetisation  for 
maximum  deformation. 


The  two-dimensional  modelling  undertaken  in  this  report  may  be  considered  qualitative  due 
to  the  (necessary)  approximation  that  the  ME  transducer,  and  also  the  magnets,  are  uniformly 
10  mm  deep  (instead  of  cylindrical).  Furthermore,  the  two-dimensional  model  assumes  that 
the  piezoelectric  response  is  isotropic,  and  mechanical  only  (i.e.  no  electromechanical  effects, 
the  piezoceramic  is  in  a  short-circuit  condition).  A  detailed  discussion  of  the  two-dimensional 
modelling  methodology  can  be  found  in  Appendix  A  (Table  A3  containing  relevant  material 
properties,  using  silver-epoxy  properties  given  in  [47]). 

2.2  Three  dimensional  modelling 

Due  to  the  cylindrical  geometry  of  the  magnets  a  three-dimensional  model  was  required  to 
obtain  reasonably  accurate  predictions  of  bearing  restoring  force  (and  also  of  strain  in  the 
piezoelectric  layer  of  the  ME  transducer,  which  is  to  be  reported  in  future  work).  Specifically, 
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the  two-dimensional  arrangement  shown  in  Figure  4  was  extended  into  three-dimensions,  as 
shown  in  Figure  5.  A  three  dimensional  model  was  developed  using  the  Comsol  application 
mode  ‘magnetostatics,  no  currents',  the  geometry  of  which  is  shown  in  Figure  5.  The  model  was 
analysed  using  Comsol  driven  by  a  Matlab  script  [48].  To  find  the  total  magnetic  force  acting 
on  the  bearing,  the  Maxwell  stress  tensor  [49]  must  be  calculated.  Comsol  employs  the 
following  projection  of  the  Maxwell  stress  tensor  onto  the  outside  surface  of  the  bearing  [43], 

nx  T2  =- pnx  -(jE.D  +  jH ,B)nx  +(nvE)DT  +(nvH)BT  (5) 


where  m  is  the  surface  normal  vector  pointing  out  from  the  surface  of  the  bearing,  p  is  the  air 
pressure,  D  =  eE  is  the  electric  displacement,  B  =  juH  is  the  magnetic  flux  density,  p  and  £  are 
the  permeability  and  permittivity  respectively.  Note  that  m,  E  the  electric  field,  and  H  the 
magnetic  field  intensity,  are  all  l-by-3  vectors.  Comsol  uses  equation  (5)  directly  in  the 
boundary  integral  of  the  stress  tensor  to  compute  the  total  force  on  the  bearing.  Magnetic  force 
predictions  were  made  assuming  that  the  ME  transducer  acted  as  a  3.5  mm  air  gap  between 
the  magnet  and  the  bearing  (i.e.  the  force  predictions  ignore  the  effect  of  the  ME  transducer). 
The  assumption  that  the  ME  transducer  acted  as  a  gap  between  the  magnet  and  the  bearing 
was  considered  reasonable  since  it  was  found  that  including  the  ME  transducer  in  the  three- 
dimensional  model  did  not  significantly  alter  the  predicted  restoring  force,  Fy.  Stable  finite 
element  solutions  were  obtained  if  the  maximum  element  dimension  chosen  was  1.5  mm  or 
less,  and  also  if  the  meshing  of  the  magnet  and  bearing  was  sufficiently  refined.  A  detailed 
discussion  of  the  three-dimensional  modelling  methodology  can  be  found  in  Appendix  B. 


Figure  5.  Three-dimensional  model  showing  the  bearing  elevated  3.5  mm  above  the  magnet  surface 


3.  Experimental 

A  simple  device  was  built  to  demonstrate  the  principle  of  bi-axial  ME  vibration  energy 
harvesting.  Figure  6a  is  a  photograph  of  the  permanent-magnet/ bearing  oscillator 
arrangement,  and  shows  the  ME  transducer  located  between  the  bearing  and  the  magnet.  The 
steel  bearing  (grade  AISI  52100)  is  spherical  with  a  diameter  of  25.4  mm;  the  rare-earth 
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magnets  (NdFeB,  grade  N42)  and  the  ME  transducer  are  cylindrical,  and  have  a  diameter  of 
10  mm. 

The  ME  transducer  was  formed  by  bonding  the  Terfenol-D  (Etrema  Products  Inc.)  and  PZT 
elements  (Pz27,  Ferroperm  Piezoceramics  A/S),  both  10  mm  diameter,  with  silver  loaded 
epoxy  (CW2400,  Circuitworks)  that  was  cured  for  two  hours  at  55  °C  under  a  compressive 
load  of  ~1  MPa.  Reference  [25]  indicates  that  the  output  from  an  ME  lamina  composite 
transducer  may  be  optimised  by  choosing  magnetostrictive  and  piezoelectric  layers  with  a 
thickness  ratio  of  approximately  three  to  one.  Hence  the  Terfenol-D  and  PZT  layers  of  the  ME 
transducer  were  chosen  to  be  1.5  mm  and  0.5  mm  thick  respectively.  The  host  structure 
oscillates  in  the  x-y  plane  (Figure  5).  The  magnets  and  the  piezoceramic  are  poled  in  the  z- 
direction.  The  magnets  are  10  mm  thick,  and  for  experimental  convenience  two  magnets  were 
stacked  together  (for  a  total  magnet  thickness  of  20  mm),  and  attached  to  a  5  mm  thick  steel 
base  using  magnetic  force.  A  thin  layer  of  beeswax  [50]  was  employed  at  the  interface  between 
the  ME  transducer  and  the  upper  magnet  to  minimise  shear  movement. 

Measurements  of  the  PZT  element  impedance  were  made  before  and  after  bonding  into  the 
lamina  composite,  using  the  approach  described  in  reference  [10].  The  results  of  the 
impedance  measurement  are  presented  in  Section  4,  however  it  is  noted  that  the  measured 
capacitance  of  ME  transducer  (i.e.  the  bonded  PZT  element)  was  1.86  nF  at  10  Hz. 

Figures  6a  and  6b  depict  the  experimental  arrangement,  showing  the  harvester  attached  to  a 
5  g  steel  host  mass  (hanging  from  two  wire  ligaments),  connected  to  a  75  N  vibration  shaker 
(TIRA  S  511  75  N).  Host  accelerations  were  measured  using  an  accelerometer  (PCB  333B40), 
and  a  laser  displacement  sensor  was  used  to  measure  the  bearing  displacement  (y-direction) 
during  testing  (ODS  115  Black-line).  Care  was  required  when  aligning  the  laser;  for  reliable 
displacement  measurements  the  laser  had  to  be  aimed  at  the  centre  of  the  bearing  parallel 
with  the  direction  of  the  host  acceleration.  A  maximum  rms  host  acceleration  of  61  mG  was 
chosen  because,  at  resonance,  larger  host  accelerations  would  drive  the  bearing  off  the  edge  of 
the  ME  transducer.  It  is  worth  noting,  as  indicated  in  Figure  3,  that  future  work  will  use  a 
cylindrical  housing  to  constrain  the  bearing  displacement  and  to  also  introduce  a  two- 
dimensional  vibro-impact  effect  that  may  lead  to  increased  operational  frequency  bandwidth 
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Figure  6.  (a)  Photo  (originally  published  in  reference  [40]),  and  (b)  schematic,  of  the  experimental 
arrangement 


UNCLASSIFIED 


9 


UNCLASSIFIED 

DSTO-TR-2649 

Two  electrical  connections  were  made  to  the  ME  transducer,  one  at  the  top  Terfenol-D 
electrode  and  one  at  the  bottom.  The  bottom  electrode  was  stainless  steel  shim  (302  grade), 
10  m  wide  and  51  pm  thick  (SHM-302-0.051x203x300,  Small  Parts  &  Bearings  Inc.)  and  was 
attached  to  the  lower  face  of  the  bottom  Terfenol-D  element  (i.e.  Figure  4)  using  double  sided 
copper  tape  ~100  pm  thick.  Wire  with  a  diameter  of  0.57  mm  was  used  for  the  top  electrode, 
and  was  attached  to  the  side  of  the  top  Terfenol-D  element  using  silver-loaded  epoxy.  After 
the  epoxy  was  cured  it  was  coated  with  silver  ink  (CW2200MTP,  Circuitworks),  this  ensured 
that  the  top  electrode  made  a  low  resistance  ohmic-contact  with  the  Terfenol-D.  Care  was 
taken  to  ensure  that  the  epoxy  did  not  interfere  with  the  motion  of  the  bearing  on  the  upper 
surface  of  the  top  Terfenol-D  element. 


(a)  (b) 

Figure  7.  Experimental  arrangements  for  the  electrical  measurements:  (a)  the  circuit  for  measuring 
open  circuit  voltage,  and  (b)  the  circuit  for  sweeping  load  resistor  R,  and  for  load  power 
measurements.  (Originally  published  in  reference  [40]). 

Two  types  of  electrical  measurements  were  made:  (i)  the  ME  transducer's  open  circuit  voltage 
(Figure  7a),  or  (ii)  the  voltage  across  a  linear  resistive  load  that  was  attached  across  the 
transducer  (Figure  7b).  Voltage  measurements  were  made  using  a  10:1  differential  probe  with 
8  MQ  input  impedance  (Hameg  HZ109),  connected  to  a  personal-computer  based  oscilloscope 
with  one  giga-sample  memory  (Picoscope  6403).  Ring-down  measurements  of  the  harvester 
(discussed  further  in  Section  4)  were  undertaken  and  suggested  that  the  resonant  frequency  of 
the  device  was  near  10  Hz.  The  measured  capacitance  of  the  ME  transducer  was  1.86  nF  (see 
above),  resulting  in  a  capacitive  reactance  at  10  Hz  of, 

Xc  =  1/(2  7r/C)~  8.6  MQ.  (6) 

This  reactance  is  larger  than  the  input  impedance  of  the  differential  probe  that  was  used  to 
measure  the  transducer  voltage,  and  hence  would  have  resulted  in  erroneous  measured 
voltages  (due  to  the  voltage  load-down  effect  of  impedances  in  parallel).  To  obtain  more 
accurate  measured  voltages  it  was  decided  to  add  90  MQ  (i.e.  nine  10  MQ  resistors)  in  series 
with  the  differential  probe,  as  shown  in  Figure  7.  To  make  allowance  for  the  voltage  divider 
effect  (of  adding  the  90  MQ  in  series  with  the  probe)  the  ME  transducer  voltages  that  were 
subsequently  measured  were  multiplied  by  d  =  (90  MQ  +  8  MQ) / (8  MQ)  =  12.25.  It  was  also 
determined  that  the  differential  probe  was  not  correctly  calibrated  at  low  frequencies  (i.e.  near 
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DC).  A 10  Hz  output  from  a  signal  generator  (50  Q  output  impedance),  set  to  a  peak  voltage 
of  2.04  V,  was  applied  to  the  differential  probe  (set  at  10:1).  The  probe  subsequently  measured 
a  peak  voltage  of  1.64  V,  indicating  that  a  further  voltage  multiplication  was  required,  being 
e  =  2.04/1.64  ~  1.244.  The  total  voltage  multiplier  was  therefore  d  x  e  =  15.24  (i.e.  voltages 
measured  by  the  differential  probe  were  multiplied  by  this  amount). 


4.  Results  and  Discussion 


The  discussion  is  divided  into  two  main  parts.  The  first  part  will  present  model  predictions  of 
the  mechanical  and  magnetostatic  behaviour  of  the  harvester.  Scaling  laws  will  also  be 
presented  that  allow  predictions  of  the  mechanical  behaviour  of  similarly  arranged  harvesters 
with  different  geometric  ratios.  The  second  part  will  discuss  the  experimental  measurements, 
both  mechanical  and  electrical,  that  have  been  made  to  determine  the  performance  of  the  bi¬ 
axial,  magnetoelectric  vibration  energy  harvester.  In  particular,  a  simple  demonstrator  has 
been  built  (see  Figure  6a)  and  its  behaviour  elucidated. 

4.1  Mechanical  and  magnetostatic  predictions 

4.1.1  Restoring  force  predictions 

Using  the  three-dimensional  model  described  in  Section  2.2  (and  Appendix  B)  force 
predictions  were  made  as  the  bearing  was  swept  in  the  y  -  z  plane  (see  Figure  5),  with  the 
origin  located  at  the  centre  of  the  uppermost  horizontal  surface  of  the  magnets. 

Varying  the  bearing  offset  in  the  z-direction  (i.e.  the  z-offset)  can  be  thought  of  as  varying  the 
thickness  of  the  ME  transducer  located  between  the  bearing  and  the  magnets.  Figure  8a  shows 
the  variation  of  the  predicted  restoring  force  Fy  acting  on  the  bearing  (in  the  y-direction)  with 
increasing  bearing  offset  in  the  ^-direction  for  a  constant  z-offset.  The  restoring  force  Fy 
appears  to  be  hookean  in  nature  about  the  centre  of  the  magnet.  The  spring  constant  (i.e. 
slope)  varies  from  1220  N/m  for  a  z-offset  of  0.5  mm  (i.e.  corresponding  to  a  transducer 
0.5  mm  thick),  to  260  N/m  for  a  z-offset  of  3.5  mm  (i.e.  corresponding  to  a  transducer  3.5  mm 
thick).  As  the  bearing  approaches  the  edge  of  the  magnet  the  restoring  force  |  Fy  |  behaves 
like  a  softening-spring,  as  shown  in  Figure  8a. 

Predictions  of  the  magnetic  attractive  force  |  Fz  |  between  the  bearing  and  the  magnet  (in  the 
z-direction)  are  shown  in  Figure  8b  as  a  function  of  bearing  offset  in  the  y-direction  (i.e.  y- 
offset).  As  might  be  expected,  the  force  increases  rapidly  as  the  bearing  moves  closer  to  the 
magnet  (i.e.  decreasing  z-offset).  For  a  constant  z-offset,  as  the  y-offset  increases  the  force 
|  Fz  |  decreases  due  to  the  cylindrical  geometry  of  the  magnet. 


UNCLASSIFIED 


11 


DSTO-TR-2649 


UNCLASSIFIED 


Figure  8.  Force  predictions  showing  (a)  \  Fy\,and(b)  \Fz\  as  a  function  of  bearing  offset  inboth  the 
y  and  z  directions 

For  the  demonstrated  harvester  (Figure  6a)  the  magnetic  force  predictions  were  made 
assuming  that  the  3.5  mm  ME  transducer  acted  as  a  3.5  mm  air  gap  between  the  magnet  and 
the  bearing.  That  case  is  represented  by  the  black  curves  in  Figure  8.  Figure  8a  predicts  that  Fy 
is  linear  near  the  central-line  (i.e.  y  <  2.0  mm)  and  has  an  effective  spring  constant  oik  = 
260  N/  m.  The  bearing  mass  is  approximately  rn  =  67  grams,  and  hence  the  predicted  resonant 
frequency  of  the  bearing/ magnet  arrangement  is, 
f RESONANCE  =  (1/  27t)n  k  /  m  ~  9.9  Hz. 
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4.1.2  Magnetic  field  predictions 

Comsol  was  also  used  to  visualise  the  static  magnetic  field  distribution  for  various  bearing 
positions.  Of  particular  interest  is  the  distribution  of  magnetic  flux  density  B  in  the  ME 
transducer  since  this  determines  the  stress  distribution  in  the  transducer,  and  hence  the 
voltage  generated.  As  described  in  Section  2.1  (and  Appendix  A)  two-dimensional  modelling 
of  the  Terfenol-D  nonlinear  magnetostrictive  response  was  carried  out  using  the  approach 
described  in  reference  [44].  Predicted  magnetic  flux  distributions  are  shown  in  Figure  9,  which 
magnifies  the  contact  region  between  the  bearing  and  the  ME  transducer.  In  Figure  9a  the 
bearing  is  centred  on  the  upper  Terfenol-D  element.  When  the  ME  transducer  is  first 
positioned  at  the  top  of  the  magnet  a  magnetostrictive  'shrinking'  effect  similar  to  that  shown 
in  Figure  2a  is  expected  to  occur,  dictated  by  the  arrangement  of  the  streamlines.  In  Figure  9b 
the  bearing  is  near  the  edge  of  the  magnet  (i.e.  4.5  mm  from  the  central-line  in  the  negative  in¬ 
direction).  The  magnetic  flux  distribution  in  the  upper  Terfenol-D  element  changes 
significantly  as  the  bearing  changes  position,  as  indicated  by  the  change  in  the  magnetic  flux 
lines  and  the  magnetic  flux  density  colour  contours  (where  red  indicates  high  density  and 
blue  low  density)  in  the  upper  Terfenol-D  element.  Again,  it  should  be  emphasised  that  the 
two-dimensional  modelling  should  be  considered  qualitative  due  to  the  (necessary) 
approximation  that  the  ME  transducer,  and  also  magnets,  are  uniformly  10  mm  deep  (i.e. 
cuboidal  instead  of  cylindrical). 
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Figure  9.  Predicted  magnetic  flux  density  and  streamlines  for  (a)  the  bearing  at  the  central  rest 
position,  and  (b)  the  bearing  near  the  edge  of  the  ME  transducer 
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Figure  10.  (a)  Predicted  magnetic  flux  density  B  in  the  y-direction  (i.e.  B.y)  as  a  function  of  distance 
along  the  lower  edge  of  the  upper  magnetostrictive  element.  B.y  is  plotted  as  a  function  of 
bearing  offset  from  the  central-line,  in  the  y-direction.  The  bearing  offset  is  varied  from  zero 
(Figure  9a)  to  a  maximum  offset  of  4.5  mm  (Figure  9b).  (b)  Mean  B.y  as  a  function  of 
bearing  offset. 

A  closer  examination  of  the  predictions  of  the  magnetic  flux  density  B,  in  the  y-direction,  for 
the  upper  Terfenol-D  element  is  presented  in  Figure  10.  As  the  bearing  moves  across  the 
surface  of  the  upper  Terfenol-D  element  (from  the  centre  to  the  edge)  Figure  10a  shows  that 
B.y  changes  direction.  The  bearing  is  steering  the  magnetic  field  through  the  ME  transducer. 
Figure  10b  presents  the  average  of  each  B.y  line  plotted  in  Figure  10a.  The  large  changes  in 
both  the  magnitude  and  direction  of  the  magnetic  flux  suggest  that  significant  strain 
variations  are  being  experienced  at  the  piezoceramic  layer. 
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Surface:  Total  displacement  [m]  Deformation:  Displacement  Max:  7.00e-7 


xl0-3  Min:  0 


Surface:  Total  displacement  [m]  Deformation:  Displacement  Max:  7.00e-7 


xlO'3  Min:  0 


Figure  11.  Predicted  deformation  of  the  ME  transducer  for  (a)  the  bearing  at  the  central  rest -position, 
and  (b)  the  bearing  near  the  edge  of  the  ME  transducer.  Deformation  has  been  amplified  by 
a  factor  of 500. 
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As  the  bearing  moves,  the  magnetic  flux  changes  (shown  in  Figures  9  and  10)  will  induce 
deformation  of  the  ME  transducer  via  the  magnetostrictive  effect.  Figure  11  shows  the 
predicted  two-dimensional  deformation  of  the  ME  transducer,  as  the  bearing  location  changes 
from  the  central  rest  position  to  near  the  edge  of  the  transducer.  The  deformation  process 
produces  mechanical  stress  on  the  PZT  element  within  the  transducer,  which  may 
subsequently  generate  a  piezoelectric  voltage.  As  already  mentioned,  the  deformations  shown 
in  Figure  11  are  indicative  only  because  of  the  geometric  approximations  in  the  two- 
dimensional  modelling. 


4.1.3  Qualitative  stress  modelling 


Shown  in  Figure  12  are  qualitative  predictions  of  the  two-dimensional  mechanical  stress  (y- 
direction)  field  in  the  ME  transducer  as  the  bearing  location  varies  from  the  central  rest 
position  (Figure  12a)  to  near  the  edge  of  the  ME  transducer  (Figure  12b).  The  predictions  are 
qualitative  only  since  they  do  not  include  any  piezoelectric  electromechanical  effects  and  also 
assume  a  cuboid  transducer  geometry,  however  they  may  provide  useful  guidance  for  future 
three-dimensional  predictions  of  stress  distributions  in  the  ME  transducer.  The  measured 
sinusoidal  open  circuit  piezoelectric  voltage  from  the  harvester  (discussed  in  detail  in 
Section  4.2)  allows  an  estimate  to  be  made  of  the  mechanical  stress  in  the  PZT  layer  of  the  ME 
transducer.  Knowing  that  the  Pz27  disk  used  in  the  ME  transducer  had  a  thickness  fp  = 
0.5  mm,  with  a  piezoelectric  voltage  constant  gsi=  -0.0107  Vm/N  [51],  and  that  the  peak  open 
circuit  voltage  was  Pout  ~  45  V  (see  Section  4.2,  Figure  15),  then  by  re-arranging  equation  (1) 
the  stress  difference  in  the  PZT  layer  can  be  estimated  as. 


V  4S 

A  erf,  =|  OUT  \= - — - T=  4.2  MPa. 

p  2xgnxtp  2x0.0107x0.5  xlO'3 


(6) 


Note  that  equation  (6)  relates  to  the  stress  variation  in  the  PZT  layer  (in  the  y-direction),  as  the 
bearing  moves  from  the  central  rest  position  to  near  the  edge  of  the  ME  transducer  (i.e. 
Figures  12a  and  12b).  However,  the  magnetic  stream  lines  depicted  in  Figure  9  suggest  that 
there  is  no  bearing  position  (on  the  surface  of  the  upper  Terfenol-D  element)  that  will  produce 
a  mechanical  stress  distribution  in  the  ME  transducer  that  is  purely  aligned  with  the  y- 
direction,  so  the  stress  difference  estimated  by  equation  (6)  will  only  be  approximate.  The 
stress  difference  estimated  using  equation  (6)  can  be  compared  with  the  predicted  stresses 
shown  in  Figure  12.  Given  that  the  average  tensile  stress  in  the  piezoelectric  element  can  be 
found  by  integrating  the  stress  in  the  layer,  then  from  Figure  12a  with  the  bearing  located  at 
the  central  rest  position,  the  average  stress  (y-direction)  in  the  piezoelectric  element  is, 

CTREST  =  6.15  N/(10*10-3*0.5*10-3)  =  1.23  MPa. 

From  Figure  12b  with  the  bearing  located  at  the  edge,  the  average  tensile  stress  (y-direction)  in 
the  piezoelectric  layer  was, 

c^edge  =  16-15  N/ (10*10-3*0.5*10-3)  =  3.23  MPa. 

The  stress  difference  in  the  piezoelectric  layer  as  the  bearing  moves  from  the  centre  to  the 
edge  is  then  ctrest-  Oedge  =  2.0  MPa,  or  about  half  the  stress  estimated  using  equation  (6)  which 
is  reasonable  considering  the  qualitative  nature  of  the  calculation.  Note  that  the  predicted 
stress  in  the  piezoelectric  layer,  in  the  z-direction,  was  negligible  for  both  bearing  locations. 
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Surface:  sy  normal  stress  global  sys.  [Pa]  Deformation:  Displacement 


Max:  7.50e6 


(b) 


Surface:  sy  normal  stress  global  sys.  [Pa]  Deformation:  Displacement 


Max:  7.50e6 


xjq-3  Min:  -5.00e6 


Figure  12.  Qualitative  two-dimensional  stress  (y-direction)  in  the  ME  transducer  for  (a)  the  bearing  at 
the  central  rest  position,  and  (b)  the  bearing  near  the  edge  of  the  ME  transducer  [i.e. 
4.5  mm  from  the  central-line] 
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4.1.4  Scaling  laws  for  the  bi-axial  vibration  energy  harvester 

For  the  harvester  arrangement  shown  in  Figure  4,  the  rest  position  for  the  bearing  is  directly 
above  the  centre  of  the  magnet  (i.e.  along  the  central-line).  This  is  not  always  the  case  since  it 
has  been  observed  that  a  small  diameter  bearing  sitting  on  the  face  of  a  large  diameter 
cylindrical  magnet  may  have  a  rest  position  away  from  the  central-line  of  the  magnet.  From 
observation,  it  is  hypothesised  that  for  generic  magnet  and  bearing  diameters,  Dmagnet  and 
Dbearing,  the  bearing  rest  position  will  remain  at  the  central-line  as  long  as  Dmagnet  <  Dbearing- 
For  a  given  ratio  D magnet/D bearing  (subject  to  the  proviso  that  Dmagnet  <  Dbearing )  the  scaling 
laws  that  apply  are: 

(i)  Frestoring  a  (Dbearing)2, 

(ii)  k  a  Dbearing, 

(iii)  /resonance  a  1/  Dbearing. 

Furthermore,  for  a  lamina  composite  of  thickness  t  (and  all  else  being  equal), 

(iv)  k  a  1/t. 

4.2  Electro-mechanical  measurements 

4.2.1  Impedance  as  a  function  of  frequency 

The  impedance  of  a  PZT  disk  with  diameter  10  mm  and  thickness  0.5  mm,  was  measured  as  a 
function  of  frequency  both  before  and  after  bonding  into  a  lamina  composite  arrangement  (as 
described  in  Section  3).  Figure  13a  shows  an  extended  frequency  range,  and  Figure  13b 
magnifies  the  region  near  the  resonant  peaks  A  and  B.  As  seen  previously  [11]  after  bonding 
of  the  PZT  disk  there  is  a  substantial  reduction  of  the  resonant  peaks  compared  with  those  of  a 
free  disk.  Equation  (6)  and  the  impedances  shown  in  Figure  13  can  be  used  to  estimate  the 
capacitance  of  the  PZT  disk;  prior  to  bonding  its  capacitance  was  2.26  nF,  and  after  bonding 
1.86  nF. 

To  confirm  that  the  90  MQ  resistor  chain  described  in  Section  3  (and  depicted  in  Figure  7)  was 
suitable  for  use  in  measuring  the  voltage  from  the  ME  transducer,  its  impedance  was 
measured  as  a  function  of  frequency.  Figure  14  shows  that  the  90  MQ  resistor  chain's 
impedance  is  constant  from  1  Hz  to  110  Hz  (with  a  small  glitch  at  the  50  Hz  mains  frequency), 
after  which  the  impedance  begins  to  roll-off,  presumably  due  to  electrical  self-resonance.  Since 
the  harvester's  mechanical  resonance  occurs  at  10  Hz,  the  110  Hz  roll-off  was  deemed 
sufficiently  high  so  as  to  not  interfere  with  the  measurement  of  transducer  voltage. 
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Figure  13.  Measured  impedance  magnitude  of  a  free  PZT  disk  with  diameter  10  mm  and  thickness 
0.5  mm,  and  of  the  disk  when  bonded  into  a  ME  transducer  arrangement  as  described  in 
Section  3.  For  the  free  disk,  the  fundamental  lateral  and  thickness  resonant  frequencies  are 
labelled  A  and  B  respectively.  Shown  in  (a)  is  the  measured  impedance  over  an  extended 
frequency  range,  and  in  (b)  a  magnification  of  the  frequency  range  near  the  resonant  peaks. 
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Figure  14.  Measured  impedance  of  the  90  MQ  resistor  chain  (shown  in  Figure  7) 


4.2.2  Magnetoelectric  voltage  and  bearing  displacements  as  a  function  of 
frequency  and  drive  angle 

The  bi-axial  nature  of  the  energy  harvester  is  clearly  seen  in  Figures  15a  and  15b.  Figure  15a  is 
a  graph  showing  nineteen  measured  frequency  sweeps  of  the  peak  open  circuit  voltage 
generated  by  the  harvester.  Sweeps  were  measured  from  low  to  high  frequency,  in  0.2  Hz 
steps.  A  maximum  rms  host  acceleration  of  61  mG  was  chosen  because,  at  resonance,  larger 
host  accelerations  would  drive  the  bearing  off  the  ME  transducer.  A  laser  displacement  sensor 
was  used  to  measure  the  bearing  displacement  (y-direction)  during  testing.  The  host 
acceleration  was  uni-axial  in  the  y-direction.  After  each  sweep  the  harvester  was  rotated  by 
5  degrees  (about  the  harvester's  central-line)  from  the  initial  position,  designated  as  angle  0°, 
through  to  an  angle  of  90°.  For  each  sweep  the  bearing's  resonant  response  occurred  at  9.8  Hz, 
similar  to  predicted  value  stated  previously.  Figure  15b  shows  the  maximum  measured  peak 
open  circuit  voltage  and  displacement  as  a  function  of  angle  (for  each  of  the  sweeps  shown  in 
Figure  15a),  and  it  can  be  seen  that  the  harvester  generates  ~  9  V  per  millimeter  of  bearing 
displacement. 
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Figure  15.  (a)  Harvester  peak  open  circuit  voltage  measured  as  a  function  of  drive  frequency,  and  (b) 
the  measured  peak  voltage  and  displacement  of  the  bearing  (at  resonance)  as  a  function  of 
drive  angle.  (Originally  published  in  reference  [40]). 
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4.2.3  Magnetoelectric  voltage  as  a  function  of  host  acceleration 

The  peak  open  circuit  voltage  of  the  harvester,  and  the  peak  displacement  of  the  bearing,  were 
measured  (at  resonance)  as  a  function  of  rms  host  acceleration.  As  the  rms  host  acceleration 
was  varied,  the  resonant  frequency  did  not  significantly  change.  Figure  16  shows  that  below 
an  rms  host  acceleration  of  40  mG  the  measured  open  circuit  voltages  (and  bearing 
displacements)  increased  in  an  approximately  linear  fashion.  There  is  however  a  clear 
gradient  increase  for  both  displacement  and  voltage  when  the  rms  host  acceleration  is  larger 
than  ~40  mG.  It  is  surmised  that  the  behavior  shown  in  Figure  16  is  due  to  the  softening¬ 
spring  effect  (as  described  previously.  Figure  8a).  It  is  worth  noting  that  the  harvester 
generates  ~  9  V  per  millimeter  of  bearing  displacement,  independent  of  the  absolute  value  of 
the  bearing  displacement. 
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Figure  16.  Measured  peak  open  circuit  voltage,  and  bearing  peak  displacement  (at  resonance),  as  a 
function  of  rms  host  acceleration.  (Originally  published  in  reference  [40]). 

4.2.4  Harvester  output  power  as  a  function  of  load  resistance 

A  sweep  through  a  range  of  load  resistor  values  was  carried  out  to  determine  the  rms  output 
power  of  the  harvester.  Using  an  rms  host  acceleration  of  61  mG  at  9.8  Hz,  the  measured  peak 
load  voltage  was  23.9  V  across  a  3.3  MQ  load  resistor  yielding  a  maximum  rms  load  power  of 
121  pW  (Figure  17).  At  the  maximum  load  power,  the  measured  peak  bearing  displacement 
was  4.61  mm,  hence  the  peak  kinetic  energy  of  the  67  gram  bearing  is  estimated  to  be  2.7  mj. 
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Figure  17.  Measured  maximum  peak  voltage,  and  rms  power,  as  a  function  of  load  resistance. 

(Originally  published  in  reference  [40]). 

4.3  Future  work 

The  modelling  of  the  bi-axial  energy  harvester,  and  laboratory  demonstration  of  the 
harvesting  concept,  both  described  in  this  report  have  provided  valuable  lessons  to  assist  with 
the  development  of  the  next  generation  of  harvester.  An  artist's  impression  of  and  assembled 
harvester  was  shown  in  Figure  3.  The  effects  of  adding  a  cylindrical  'bearing  housing'  to 
constrain  the  bearing's  movement  and  also  to  introduce  a  vibro-impact  effect  [9]  will  need  to 
be  examined.  Rotational  effects  produced  by  the  bearing  rolling  across  the  surface  of  the 
Terfenol-D  were  ignored  in  this  report,  and  will  need  to  be  examined  in  future  work.  For 
example  the  rotational  effects  may  prove  important  in  addressing  surface  wear  at  the 
bearing/Terfenol-D  interface  which  has  been  found  to  be  an  issue.  The  modelling  presented 
in  this  report  is  preliminary.  Of  particular  importance  is  the  inclusion  of  piezoelectric  (i.e. 
electromechanical)  effects  in  the  models,  which  will  require  the  use  of  the  Comsol 
piezoelectric  application  mode.  This  will  allow  examination  of  the  electric  fields  in  the 
piezoelectric  layer,  and  hence  allow  predictions  to  be  made  of  the  magnetoelectric  voltages 
being  generated  by  the  ME  transducer.  The  modelling  of  the  Terfenol-D  layers  will  need  to  be 
more  complex  to  address  the  magnetostrictive  response  at  low  applied  magnetic  field  levels, 

in  particular  for  the  situation  /  IIS<  1/V3  in  equation  (4).  Reference  [44]  adopts  a 
magnetostriction  constant  of  As  =2  x  1CT4,  the  validity  of  this  choice  of  As  [used  in  equation 
(4)]  will  be  further  examined.  To  fully  understand  the  bi-axial  harvesters  operation  the  three- 
dimensional  stress  (and  strain)  distribution  within  the  ME  transducer  needs  to  be  examined  as 
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a  function  of  bearing  position.  Eddy  current  losses  need  to  be  quantified,  although  they  are 
expected  to  be  small  due  to  the  low  frequency  nature  of  the  bi-axial  harvester.  Finally,  the 
response  of  the  bi-axial  harvester  to  a  bi-axial  excitation  needs  to  be  examined. 


5.  Conclusion 


DSTO  is  investigating  vibration  energy  harvesting  approaches  with  the  ultimate  goal  of 
capturing  airframe  vibrational  energy,  and  then  converting  it  into  electrical  power  for  use  by 
structural  health  monitoring  sensors  and  devices.  This  report  investigated  a  bi-axial  vibration 
energy  harvesting  approach  with  the  aim  of  harvesting  energy  from  the  multi-axis 
accelerations  that  can  be  experienced  by  an  aircraft.  A  bi-axial  oscillator  was  created  using  a 
permanent-magnet/ ball-bearing  arrangement.  The  magnet  produces  a  bi-axial  restoring  force 
on  the  bearing,  and  as  the  bearing  oscillates  it  steers  magnetic  field  through  a 
magnetostrictive/  piezoelectric  laminate  transducer  thereby  producing  an  oscillating  charge 
that  can  be  harvested.  The  chosen  harvester  arrangement  consisted  of  a  steel  ball-bearing  with 
25.4  mm  diameter,  and  a  rare-earth  magnet  and  magnetoelectric  (ME)  transducer  both  with 
10  mm  diameter.  Modelling  indicated  that  the  magnetic  restoring  force  has  an  effective  spring 
constant  of  260  N/m  when  the  ME  transducer  3.5  mm  thick.  Modelling  was  also  used  to  make 
a  qualitative  assessment  of  the  magnetic  flux  changes  in  the  ME  transducer  as  the  bearing 
oscillates,  which  indicated  that  large  flux  variations  occur  as  the  bearing  moves  from  the 
magnets  central-line  towards  the  edge.  A  simple  laboratory  demonstrator  of  a  bi-axial  ME 
energy  harvester  was  created  using  a  Terfenol-D/lead  zirconate  titanate/Terfenol-D 
transducer.  Flarvester  output  was  measured  as  a  function  of  drive-angle,  host  acceleration, 
and  load  resistance.  The  harvester  produced  a  peak  rms  power  of  121  pW  from  an  rms  host 
acceleration  of  61  mG  at  9.8  Hz. 
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Appendix  A:  Two-dimensional  Modelling 


This  appendix  provides  detailed  instructions  on  using  Comsol  multi-physics  software  [42]  to 
carry  out  two-dimensional  modelling  of  the  bi-axial  magnetoelectric  vibration  energy 
harvester  described  in  this  report.  Two-dimensional  modelling  of  the  harvester's 
magnetostrictive  (Terfenol-D)  elements  follows  the  approach  outlined  in  reference  [44]. 

MODEL  NAVIGATOR 

1.  Start  COMSOL  Multiphysics  software. 

2.  Under  Space  Dimension,  choose  2D  (see  Figure  Al). 

3.  In  the  application  mode  list,  open  AC/DC  Module>Statics, 
Magnetic>Perpendicular  Induction  Currents,  Vector  Potential. 

4.  Click  the  Multiphysics  buttonand  then  the  Add  button. 

5.  In  the  application  mode  list,  open  Structural  Mechanics  Module>Plane  Strain. 

6.  Click  the  Add  button. 

7.  Click  OK  to  close  the  Model  Navigator. 


Figure  Al.  Comsol  model  navigator 
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GEOMETRY  MODELLING 

1.  Shift-click  the  Rectangle/Square  button  on  the  Draw  toolbar  for  each  shape,  and 
enter  the  parameters  given  in  Table  A1  (the  description  is  only  for  the  readers 
information). 


Table  Al. _ Geometric  settings  for  the  harvester  model 


Size 

Position 

Width 

Height 

X 

Y 

(Description) 

0.3 

0.2 

-0.15 

-0.1 

Air 

0.07 

0.005 

-0.035 

-0.029 

Base  Plate 

0.01 

0.012 

-0.005 

-0.024 

Lower  Magnet 

0.01 

0.012 

-0.005 

-0.012 

Upper  Magnet 

0.01 

1.0E-4 

-0.005 

0 

Bond 

0.01 

0.0015 

-0.005 

1.0E-4 

Lower  Terfenol-D  layer 

0.01 

2.5E-5 

-0.005 

0.0016 

Bond 

0.01 

5.0E-4 

-0.005 

0.001625 

PZT  layer 

0.01 

2.5E-5 

-0.005 

0.002125 

Bond 

0.01 

0.0015 

-0.005 

0.00215 

Upper  Terfenol-D  layer 

2.  Shift-click  the  Ellipse/Circle  (Centred)  button. 

3.  In  the  Circle  dialog  box,  type  0.0127  in  the  Radius  edit  field  and  set  the  centre 
coordinates  to  (0,  0.1635)  by  typing  0  in  the  X  edit  field  and  0.01635  in  the  Y  edit 
field  (this  is  the  ball  bearing). 

4.  Click  the  Point  button  on  the  Draw  toolbar,  and  draw  a  point  at  (0,  0).  The 
completed  geometry  should  look  similar  to  that  shown  in  Figure  A2. 
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OPTIONS 

1.  From  the  Options  menu,  open  the  Constants  dialog  box  (see  Figure  A3). 

2.  Enter  a  constant  with  the  name  Lambda_s,  the  expression  200E-6,  and  the 
description  magnetostriction  constant  (the  description  is  optional). 

3.  Enter  a  constant  with  the  name  Msat,  the  expression  1500E3[A/m],  and  the 
description  saturation  magnetisation  (the  description  is  optional). 

4.  Click  Apply  and  OK  to  close  the  Constants  dialog  box. 


Figure  A3.  Comsol  constants  dialog  box 

5.  From  the  Options  menu,  open  the  Functions  dialog  box. 

6.  Click  the  New  button. 

7.  In  the  Function  Name  field,  enter  HBFe. 

8.  Select  the  Interpolation  tick  box,  and  select  File  in  the  list  to  Use  data  from. 

9.  Locate  the  data  file  required,  or  alternatively  enter  the  data  as  a  table  (the 
multimedia  accompaniment  to  this  report  contains  the  file  'HBFe.txt',  also  see 
Figure  A4). 

10.  Select  Linear  in  the  Interpolation  method  list. 

11.  Click  Apply  and  OK  to  close  the  Functions  dialog  box. 
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Defined  functions 


New. 


Function  definition 

Function  name:  HBFe 


Interpolation  method:  Linear 
Extrapolation  method:  interpolation  function 
Value  outside  range: 

Use  space  coordinates  as  default  function  arguments 

Function  description _ 


a 

a 


r  <$  £>  &  jf  lo,  ,4 


HBFe(x) 


□ 

□ 


Figure  A4.  Comsol  functions  dialog  box 


12.  From  the  Options  menu,  under  Expressions,  select  Subdomain  Expressions. 

13.  Select  subdomains  6  and  10,  and  an  expression  with  the  name  Lambdajx,  and  the 
expression  L5*Lambda_s*(Mx_emqa/Msat)A2.  (See  Figure  A5). 

14.  Then  enter  another  expression  with  the  name  Lanibda_y,  and  the  expression 
L5*Lambda_s*(My_emqa  /  Msat)  A2. 
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Figure  A5.  Comsol  subdomain  expressions  dialog  box 

PHYSICS  SETTINGS 

Subdomain  Settings  -  Perpendicular  Induction  Currents,  Vector  Potential 

1.  From  the  Multiphysics  menu,  select  the  Conductive  Media  DC  application  mode 
(see  Figure  A6). 

2.  From  the  Physics  menu,  select  Subdomain  Settings. 

3.  In  the  Subdomain  Settings  dialog  box,  select  each  subdomain  and  enter  the 
parameters  given  in  Table  A2. 


Table  A2.  Subdomain  settings  for  the  harvester  model 


Subdomains 

Library 

Material 

L 

H<->B 

Other 

1,  5,  7,  8,  9 

NA 

0.01 

B  =  p0prH 

Pr  =  1 

2,11 

Soft  Iron 

(with  losses) 

0.01 

NA 

NA 

3,4 

NA 

0.01 

B  =  p0prH  +  Br 

Pr  =  1 

B,-  =  (0, 1.3) 

6, 10 

NA 

0.01 

H  =  f ( |  B  |  )eB 

|  H  I  = 

HBFe(normB emqa[l/T])[A/m] 

4.  Click  Apply  and  OK  to  close  t 

re  Subdomain  Settings  dialog  box. 
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Figure  A6.  Comsol  subdomain  settings  dialog  box 

Subdomain  Settings  -  Plane  Strain 

1.  From  the  Multiphysics  menu,  select  the  Plane  Strain  application  mode  (see  Figure 
A  7). 

2.  From  the  Physics  menu,  select  Subdomain  Settings. 

3.  In  the  Subdomain  Settings  dialog  box,  select  subdomains  1,  2,  3, 4, 11  and  deselect 
the  tick  box  labeled  "Active  in  this  domain". 

4.  For  the  other  subdomains  enter  the  parameters  given  in  Table  A3. 


Table  A3.  Material  properties  and  geometric  settings  for  the  harvester  model 


Subdomains 

Library  Material 

E 

V 

P 

Thickness 

5 

NA 

50E6 

0.33 

958 

0.01 

7,9 

NA  (see  ref.  [47]) 

2.9E6 

0.35 

3890.2 

0.01 

8 

Lead  Zirconate 

Titanate  (PZT-5A) 

59E9 

0.33 

NA 

0.01 

6, 10 

NA 

30E9 

0.45 

7850 

0.01 

5.  For  subdomains  6  and  10,  click  the  Initial  Stress  and  Strain  button  and  select  the 


tick  box  labeled  "Include  initial  strain". 

6.  For  Initial  normal  strain  type  ( Lambda_x ,  Lambda_y,  0). 

7.  Click  Apply  and  OK  to  close  the  Subdomain  Settings  dialog  box. 
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Figure  A7.  Comsol  subdomain  settings  dialog  box,  initial  stress  and  strain  menu  tab 


Boundary  Settings  -  Plane  Strain 

1.  From  the  Physics  menu,  select  Boundary  Settings  (see  Figure  A8). 

2.  In  the  Boundary  Settings  dialog  box,  select  boundaries  12  and  24  and  select 
Prescribed  Displacement  in  the  Constraint  Condition  list. 

3.  Under  Standard  notation,  select  the  Ry  tick  box,  and  verify  the  field  displays  0. 

4.  Click  Apply  and  OK  to  close  the  Boundary  Settings  dialog  box. 
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Figure  A8.  Comsol  boundary  settings  dialog  box 

Point  Settings  -  Plane  Strain 

1.  From  the  Physics  menu,  select  the  Point  Settings  (see  Figure  A9). 

2.  In  the  Point  Settings  dialog  box,  select  point  15. 

3.  Under  Standard  Notation,  select  the  Rx  and  Ry  tick  boxes,  and  verify  that  the 
fields  display  0. 

4.  Click  Apply  and  OK  to  close  the  Point  Settings  dialog  box. 


Figure  A9.  Comsol  -point  settings  dialog  box 
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MATERIALS  PROPERTIES 

1.  From  the  Option  menu,  open  the  Materials/Coefficients  Library  dialog  box. 

2.  From  the  Materials  area,  expand  the  Model  heading  and  select  Soft  Iron  (with 
losses). 

3.  From  the  All  tab,  scroll  down  the  list  of  coefficients  and  find  sigma. 

4.  Clear  the  sigma  edit  field  and  type  0. 

5.  Click  Apply  and  OK  to  close  the  Materials/Coefficients  Library  dialog  box. 
MESH  GENERATION 

Subdomain  Properties 

1.  From  the  Mesh  menu,  select  Free  Mesh  Parameters. 

2.  In  the  Free  Mesh  Parameters  dialog  box,  select  the  Subdomain  tab. 

3.  Select  each  subdomain  and  enter  the  paramaters  given  in  Table  A4. 

Table  A4.  Mesh  settings  for  harvester  model. 


Subdomains 

Maximum  element  size 

Method 

1, 11 

NA 

Triangle 

2,3,4 

IE-3 

Triangle  (advancing  front) 

5,  6,  7,  8,  9, 10 

IE-4 

Triangle  (advancing  front) 

Boundary  Properties 

1.  In  the  Free  Mesh  Parameters  dialog  box,  select  the  Boundary  tab. 

2.  Select  boundaries  37  and  39,  and  enter  3E-4  in  the  Maximum  element  size  edit 
field. 

3.  Select  boundaries  23  and  25,  and  enter  IE-4  in  the  Maximum  element  size  edit 
field. 

4.  Click  the  Remesh  button  to  generate  the  mesh  (as  indicated  in  Figure  A10). 

5.  Click  Apply  and  OK  to  close  the  Free  Mesh  Parameters  dialog  box. 


Figure  A10.  Comsol  meshing  of  the  two-dimensional  model 
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COMPUTING  THE  SOLUTION 

1.  From  the  Solve  menu,  select  Solver  Parameters  (see  Figure  All). 

2.  In  the  Analysis  types  area,  select  Static  in  the  Perpendicular  Induction  Current, 
Vector  Potential  list. 

3.  Verify  that  in  the  Analysis  types  area  the  Plane  Strain  list  is  also  Static. 

4.  Verify  that  in  the  Linear  system  solver  area  and  list.  Direct  (UMFPACK)  is 
selected. 

5.  Click  Apply  and  OK  to  close  the  Solver  Parameters  dialog  box. 

6.  Click  the  Solve  button  (equal  sign)  on  the  Main  toolbar  to  solve  the  model  (or 
choose  Solve  Problem  from  the  Solve  menu).  Example  output  is  shown  in  Figure 
A12. 


Figure  All.  Comsol  solver  parameters  dialog  box 
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POSTPROCESSING  AND  VISUALISATION  EXAMPLE 


Figure  All.  Screen  capture  using  Comsol  processing  and  visualisation 
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Appendix  B:  Three-dimensional  modelling 

This  appendix  provides  specific  details  on  utilising  Comsol  multi-physics  software  [42]  to 
carry  out  three-dimensional  modelling  of  the  bi-axial  magnetoelectric  vibration  energy 
harvester  described  in  this  report.  Three-dimensional  modelling  of  the  harvester's 
magnetostrictive  (Terfenol-D)  elements  is  done  by  adapting  the  two-dimensional  approach 
outlined  in  reference  [44] . 

Comsol  was  used  to  model  the  three-dimensional  forces  on  the  ball-bearing,  with  the 
“Magnetostatics,  No  Currents"  application  mode  from  the  “AC/DC  Module".  Figure  B1  shows  an 
example  of  the  geometry  and  mesh.  Both  the  large  cylinder  and  the  large  sphere  are 
modelling  air  as  a  magnetic  conducting  environment,  and  the  cylinder  has  a  finer  mesh  for 
greater  accuracy. 


Figure  Bl.  Three-dimensional  geometry  and  mesh  used  for  modelling. 

The  Matlab  interface  was  used  to  automate  a  sweep  of  different  ball  bearing  positions.  The 
vertical  air  gap  represents  the  height  of  the  transducer  and  the  ball  bearing  was  swept  along 
the  y-axis  to  develop  the  restoring  force  curve  at  multiple  heights. 

The  script  for  these  automatic  sweeps  is  included,  which  repeatedly  alters  the  geometry 
accordingly,  generates  a  new  mesh  for  this  geometry,  then  solves  the  model  and  records  the 
forces.  The  ball  bearing  and  air  cylinder  surrounding  it  are  both  moved  in  each  cycle,  while  all 
other  elements  remain  fixed.  The  mesh  consisted  of  around  33000  elements  and  was  solved  for 


UNCLASSIFIED 


41 


UNCLASSIFIED 

DSTO-TR-2649 

around  146000  degrees  of  freedom;  on  a  machine  with  16  GB  RAM  and  a  2.67  GHz  CPU  a 
sweep  of  91  instances  took  about  350  minutes  to  run. 

THREE-DIMENSIONAL  MODELLING  SCRIPT 


%  JEM_2DpositionSweep 

%  Uses  3D  magnetic  forces  model  to  solve  for  a  sweep  of  ball  bearing 
%  positions  in  a  2D  plane. 

% 

%  ©version  3.1  -  17/05/2011 
%  ©author  Josh  McLeod 

zShiftSet  =  . 5:. 5:3.5; 
yShiftSet  =  0:.5:4.5; 

totalSolNo  =  length (zShiftSet)  *  length (yShiftSet ) ; 
clc; 

j  =  0; 

solNo  -  0; 

data2  =  {'zShift',  'yShift',  'Fy',  'Fz'}; 

fprintf ( ' \t\t\t--  2D  POSITION  SWEEP  — \n\t\ t\ t\ t\ t Josh  McLeod\n\n ' ) ; 

%  COMSOL  version 
clear  vrsn 

vrsn.name  =  'COMSOL  3.5'; 
vrsn . ext  =  ' a ' ; 
vrsn. major  =  0; 
vrsn. build  =  603; 
vrsn. res  =  '$Name:  $'; 

vrsn. date  =  '$Date:  2008/12/03  17:02:19  $'; 
fem. version  =  vrsn; 

gl=sphere3 ( ' . 1 ' , ' pos ',{'0', 'O', '0.0125'}, ' axis ',{'0','0','l'},' rot ' , ' 0 ' ) ; 

g3=block3 ( ' . 1 ' ,  '  .  1 '  ,  '  .  005  '  ,  'base ' ,  ' centre ' ,  'pos ', {'O',  'O',  '  —  . 0025 ' } , ' axis ',{'0','0','l'},' rot ' ,  ' 0 ' ) ; 
g4=cylinder3 ( ' . 005 ' , ' . 02 ' , 'pos ',{'0','0','0'},' axis ',{'0','0','l'},' rot ' , ' 0 ' ) ; 

%  Z  position  loop 
for  zShift  =  zShiftSet 

j  -  j  +  i; 

i  =  0 ; 

data(j+l,l)  =  {['zShift  =  ',  num2str ( zShift) ]} ; 

fprintf (' \nSolving  for  zShift  =  %f\n',  zShift); 

%  Y  position  loop 
for  yShift  =  yShiftSet 

i  =  i  +  1; 

solNo  =  solNo  +  1; 

if  (j==l) 

data(l,i+l)  =  {['yShift  =  ',  num2str (yShift) ]} ; 
if  (i==l) 

tic;  %  start  timer  to  check  first  solve  time 

end 

end 

fprintf (' \tyShift  =  %f\n',  yShift); 

%  Geometry 
clear  draw 

g2=cylinder3 ( ' . 015 ' ,  '  . 05 ' ,  'pos ' , { ' 0 ' ,  (yShift/1000) ,  ' 0 ' } , ' axis ',{'0','0','l'},  ' rot ' ,  ' 0 ' ) ; 
g5=sphere3 ( ' .0127', 'pos', {'O', (yShift/1000) , ( . 020+ . 0127+zShif t/1000 ) },  ... 

' axis ',{'0','0','l'},' rot ' ,  ' 0  ' ) ; 

%  Analyzed  geometry 
clear  s 

s .obj  s={gl, g2, g3, g4, g5} ; 
f em. draw=struct ( ' s ' ,  s)  ; 
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f  em .  geom=geomcsg  ( fem)  ; 

%  Initialize  mesh 
f em.mesh=meshinit (fem,  .  .  . 

' hauto ' ,  9,  ... 

' hmaxfac' ,  [14,1.56-3, 16, 1 .5e-3, 21,1 .5e-3, 25, 1 .5e-3,  30, 1 .5e-3]  ,  ... 

' hmaxsub ' , [1, le-2,2, le-2, 3, 1 . 5e-3, 4 , 1 . 5e-3, 5, 1 . 5e-3] ) ; 

%  Application  mode  1 
clear  appl 

appl .mode . class  =  'MagnetostaticsNoCurrents ' ; 

appl. module  =  'ACDC'; 

appl . assignsuf fix  =  '_emnc' ; 

clear  bnd 

bnd.type  =  { ' VmO ' , ' cont ' } ; 
if  (yShift==0) 

bnd.ind  =»  [1,1,1,  1,2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,1,  1,2, 2, 2, 2,  1,2, 2, 
1 , 2 , 2 , 2  ]  ; 

else 

bnd.ind  =  [1,1, 1,1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,1, 2, 2, 2, 2, 1,1, 2, 
2 , 2 , 2 , 2  ]  ; 

end 

appl. bnd  =  bnd; 
clear  equ 

equ .magconstrel  =  {'mur', 'mur', 'mur', 'Br'}; 

equ .M  =  {  { 0 ; 0 ; 0 } , { 0 ; 0 ; 0 } , { 0 ; 0 ; 0 } , { 0 ; 0 ; 10e6 } } ; 

equ. mur  =  {1,100,100,1}; 

equ .maxwell  =  {{},{},' BB ' , ' MM ' } ; 

equ . Br  =  { { 0;  0;  0 } ,  { 0;  0; 0 } , { 0; 0; 0 } , { 0; 0; 1 . 3 }  }  ; 

equ.ind  =  [1,2, 1,3, 4]; 

appl. equ  =  equ; 

fem.appl{l}  =  appl; 

fem. frame  =  {'ref'}; 

fem. border  =  1; 

clear  units; 

units .basesystem  =  'SI'; 
fem. units  =  units; 

%  ODE  Settings 
clear  ode 
clear  units; 

units .basesystem  =  'SI'; 
ode. units  =  units; 
fem.ode=ode; 

%  Multiphysics 
f em=multiphysics (fem)  ; 

%  Extend  mesh 

f em . xmesh=meshextend ( fem,  ... 

' dofversion ' ,  1 ) ; 

%  Solve  problem 
f em . sol=f emstatic (f em,  ... 

' solcomp ' , { ' Vm' } ,  ... 

' outcomp ' , { ' Vm ' } ,  ... 

' blocksize ' , ' auto ' ,  ... 

' linsolver ' , ' eg ' ,  ... 

' prefun ' , ' amg ' ) ; 

%  Solution  data 

force (1 ) =postint ( fem, ' BB_forcey_emnc ' ,  ... 

' unit ' , 'N ' ,  ... 

' recover ',' of f '  ,  ... 

'dl',23,  ... 

' edim ' , 0) ; 

force (2 ) =postint ( fem, ' BB_f orcez_emnc ' ,  ... 

' unit ' , ' N ' ,  ... 

' recover ',' of f ' ,  ... 

'dl',23,  ... 

' edim ' , 0) ; 

%  Two  methods  of  storing  data 
data ( j+1, i+1)  =  {force}; 

data2 (solNo+1 , : )  =  {zShift,  yShift,  force (1),  force(2)}; 
fprintf ( ' \t\t (Fy  =  %f  and  Fz  =  %f)\n',  force (1),  force (2)); 
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%  Calculate  expected  completion  time 
if  ( j  ==1  &&  i  ==1) 

currentSolveTime  =  toe; 

f ullSolveTime  =  currentSolveTime  *  totalSolNo; 

%  Initialise  progress  bar 

progress  =  waitbar (solNo/totalSolNo,  ... 

['Test  will  run  for  approximately  int2str ( fullSolveTime/ 60 ) ,  ... 

'  minutes  (has  been  running:  '  int2str (currentSolveTime/60) ,  '  minutes) ']); 

else 

%  Update  progress  bar 
currentSolveTime  =  toe; 

waitbar (solNo/totalSolNo, progress,  . . . 

['Test  will  run  for  approximately  int2str ( totalSolNo/solNo*currentSolveTime/ 60 ) , 

'  minutes  (has  been  running:  '  int2str ( toe/ 60 ) ,  '  minutes)']); 

end 

end  %  Y  position  loop 
end  %  Z  position  loop 

waitbar (1 , progress, [' Finished  in  int2str (toc/60)  '  minutes']); 
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